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We investigate the formation of ripples on the surface of wind- 
blown sand based on the one-dimensional model of Nishimori and 
Ouchi [Phys. Rev. Lett. 71, 197 (1993)], which contains the pro- 
cesses of saltation and grain relaxation. We carry out a nonlinear 
"■Q | analysis to determine the propagation speed of the restabilized rip- 

Ph ' pie patterns, and the amplitudes and phases of their first, second, 

and third harmonics. The agreement between the theory and our nu- 
merical simulations is excellent near the onset of the instability. We 
also determine the Eckhaus boundary, outside which the steady ripple 

t> ■ patterns are unstable. 

in 

^ ■ PACS numbers: 47.54.+r, 45.70-n, 92.10.Wa 

§^ : 1 Introduction 

-i— > 

Since the pioneering work of Bagnold [1] , many researchers have investigated 
the complex dynamics of dry granular materials at a surface [2-6] . Dry gran- 
ular materials are assemblies of macroscopic objects that interact with each 
other essentially via a hard core repulsive potential. Hence they are loosely 
connected, particularly at the surface. When those grains at the surface are 
exposed to a wind, they can readily be ejected and carried by the wind until 
gravity eventually pulls them back to the surface. The dynamics of a single 
grain is rather simple, given by the Newtonian trajectory of a point particle. 
Even so, experiments have shown that the collective response of the grains 
can become exceedingly complex, ranging from formation of simple ripple 
patterns to ridges and dunes to violent tornadoes [2]. Our current under- 
standing of such complex phenomena remains mostly confined to compiling 
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data on experimental observations. With regard to the formation of ripple 
patterns, however, there have been some attempts to construct a simple yet 
physical continuum model. 

We will investigate the continuum model due to Nishimori and Ouchi 
[3]. The Nishimori-Ouchi (NO) model of ripple patterns accounts for two 
elementary processes of sand transportation by the wind which have been 
identified by investigators in aeolian sand dynamics, namely saltation and 
creep [1,2]. Saltation refers to the process by which surface grains are ejected 
into the air under the influence of a strong wind, and are blown downwind 
where they collide with other surface grains. There they transfer momentum 
to these downwind grains, which may themselves be ejected in turn, thereby 
continuing the process (see Fig. 1). Creep is the surface movement of grains 
too heavy to be ejected into the air but light enough to be pushed along the 
surface. Creep also describes the surface movement of grains on hills under 
the influence of gravity. Previous studies based on the NO model have been 
confined largely to linear stability analysis and Monte Carlo simulations of a 
lattice version with simple rules for the grain dynamics [3,11]. The purpose 
of this paper is to go one step farther by carrying out a nonlinear analysis of 
the continuum model and uncovering some of the features of ripple formation 
that are inaccessible to linear analysis. 

In particular, we carry out a weakly nonlinear analysis valid near the 
onset of instability of a flat sandbed to determine the amplitude, shape, and 
propagation speed of the ripple pattern that forms in this regime. We also 
compare these results with our numerical integrations of the model equa- 
tions. These computations are rather unusual because the model lacks an 
up-down symmetry, and especially because accounting for saltation makes 
the model nonlocal in space. We find, however, that the process of pattern 
selection in this simple one-dimensional system, in particular the selection of 
the wavelength and speed of the patterns [7], is similar to what is seen in more 
complicated multidimensional systems such as directional solidification [8] or 
directional viscous fingering [9]. That is, the wavelength of the final pattern 
depends on the initial conditions, and may lie anywhere within a band of 
linearly stable final states. The stable band turns out to be somewhat wider 
than in most other models. 

In the next section we review the Nishimori-Ouchi model equations [3], 
point out a physical symmetry which they violate, and propose a simple 
modification of the model which respects that symmetry. In Section III we 



carry out a linear stability analysis of the flat-sandbed solution of the model, 
both for the original Nishimori-Ouchi equations and for our modification. We 
extend this in Section IV to give a weakly nonlinear analysis for both forms of 
the model. Section V presents our numerical calculations and compares them 
with the results of the weakly nonlinear analysis. The results are discussed 
in the final section. 

2 One-Dimensional Model for Windblown Sand 

The starting point of the Nishimori-Ouchi (NO) model [3] is a local conser- 
vation law for sand grains. Let h(x, i) be the local height of the sand bed 
at position x and time t, measured from some reference level. The height 
increases when grains are added at position x. We write 

m + ~dx~ ~ Qnh (1) 

where Ji(x,t) is a local flux of grains in the positive direction at x, and 
Qni(x,t) is the net input of grains at x due to nonlocal processes. 

The expression for J\ embodies a model of creep. Nishimori and Ouchi 
choose J i = —D(dh/dx). Note that this merely expresses the tendency of 
grains to roll downhill; it does not include any bias favoring motion in the 
direction of the wind. 

Saltation is modeled by gain and loss terms in the nonlocal transfer rate 
Q n i. Let N(x, t) denote the outward saltation flux of particles from x at time 
t. That is, let N(x, t) dx be the number of particles per unit time taking off 
from positions between x and x+dx. The loss term in Q n i is then —A N(x, t), 
where A is a scale parameter. The gain term is proportional to the rate at 
which particles arrive at x from other locations £ upwind of x. Suppose all 
particles which take off from the interval (£, £ + d£) subsequently land in the 
interval (x,x + dx). Then the number of particles per unit time landing in 
this latter interval of length dx is iV(£, t) d£, so the gain term in the saltation 
flux is then AN(£,t) (d^/dx). It is possible to have more than one £ which 
satisfies this equation for a given x. That is, grains landing at x may have 
come from more than one takeoff point £. If this is the case, then the input 
term in the evolution equation should be summed over the different values 
of£. 



Note that evaluating N(£) at time t neglects the flight time of the incom- 
ing grains; we expect the evolution of the sandbed profile to take place on 
a much longer time scale than this, so that the time delay between takeoff 
and landing should be unimportant. Indeed, experiments on ripple formation 
by sand transported by water [10] show the evolution of the ripple pattern 
occurring on time scales of several hours. 

Combining the various contributions to the flux and substituting into the 
general conservation law for h gives the model evolution equation for the 
sandbed profile, 



dh d dh 
dt dx dx 



N(U)^-N(x,t) 



(2) 



Note that this equation is nonlocal in x, as a result of the saltation gain 
term, which depends on conditions at a position £ which is a finite distance 
upwind of x. 

To complete the model, we must now specify the saltation function, an 
equation for the flight length of a single grain. In general, the locations x 
and £ in the evolution equation will be related by x = £ + L, where L is 
the horizontal length an ejected grain travels from takeoff to landing. This 
will depend on the size of the grain, its speed when it takes off, the wind 
velocity profile, and the topography of the sandbed itself. Nishimori and 
Ouchi proposed the simple ansatz 

L = L + bh(£,t)- (3) 

Here, Lq is a parameter proportional to the shear stress of the wind at the 
surface, or more precisely to the friction velocity of the wind on the sand 
surface [11], and h in general depends on the average drag force on the grain. 
Nishimori and Ouchi took both L and b to be constant, essentially assuming 
the wind velocity to be a constant, independent of x and t and unaffected by 
changes in the sandbed profile. 

Equation (H) merely indicates that the higher the takeoff point of a grain 
in saltation, the longer its trajectory. As Nishimori and Ouchi point out 
[3], this amounts to assuming that the height and topography at the point 
of landing may be neglected, and that only the surface height (as opposed 
to local topography) is important at the takeoff point. While this may be 
reasonable if h(x) is everywhere close to zero, it does violate a symmetry 



of the physical problem, namely that the dynamics should be unaffected if 
we add any constant to h, thus changing our reference level. To restore this 
symmetry, it may be more appropriate to take the saltation function to be 

L = Lo + b[h(£)-h(x)], (4) 

where £ is the takeoff point and x is the landing point. We will discuss the 
effects of this modification below. 

For convenience, we now put the model into dimensionless form. Taking 
L , b and D to be constants, we choose L to be the unit of horizontal 
length, L /b to be the unit of vertical length (i.e., of h), and L\jD to be the 
time unit. Further, we define J(x,t) = (AbL Q /D) N(x,t), a dimensionless 
measure of the outward grain flux due to saltation. With these definitions, 
the evolution equation (§) becomes 

dh d 2 h df 

with the original NO saltation relation becoming the condition 

x = Z + l + h(£,t)- (6) 

The model simplifies further if we choose J(x, t) to be a constant J inde- 
pendent of x and t, an assumption whose physical content is that the wind 
is uniform and there is no flux dependence on surface height. The evolution 
equation is then 

at dx 2+ \dx )■ [ ] 

This is the form of the problem which we will analyze below, using both the 
NO saltation relation @ and our symmetric modification of it, 

x = £+l+h{Z,t)-h{x,t). (8) 

3 Linear Stability Analysis 

We first note that a flat sandbed, h = ho = constant, is always a steady- 
state solution of the model, for either choice of saltation relation. For the 
symmetric saltation relation (M) this always gives £ = x — 1, while for the 



NO relation (|6|) we have £ = x — 1 — /to- In the latter case, however, we may 
then redefine the length and time units - and the value of J - to map the 
solution with any finite h Q (provided h > — 1) onto the solution with h Q = 0. 
Specifically, we would take the horizontal length unit to be L (l + h ) instead 
of L , and J would then be AbL Q (l + ho)N/D rather than AbL N/D. Thus 
we will take h = to be the steady state whose stability we will investigate. 
When h is small, we may linearize the NO saltation relation to get 



£ ~ x — 1 — h(x — 1), 



(9) 



so that there is a single, unique £ for each x. From this we obtain d^/dx = 
1 — h'(x — 1), where the prime indicates partial differentiation with respect 
to x. The linearized evolution equation is then 



dh 

^ = h"{x,t)-Jh'{x-l,t). 



(10) 



Linear stability analysis proceeds in the usual way: we write h(x, t) as a 
linear combination of Fourier modes, 



h(x,t) = / hk(t) exp(ikx)dk, 



(11) 



substitute this into the linearized evolution equation, and note that - even 
with the nonlocal term present - the modes do not couple. Thus we find that 
each mode grows or decays exponentially with time, 



with [3,11] 



h k (t) ocexp[(cr fc -iu k )t], 



Ck = —k — J k sink 
u k = J k cos k. 



(12) 



(13) 
(14) 



If we use the symmetric saltation relation instead of the NO relation, then 
the linearized evolution equation has an extra term +Jh'(x,t) on the right 
side. This leaves a k unchanged, but replaces ([14]) by ui k = — J k (1 — cos A;). 
The locus a k = in the J-k plane defines a stability boundary. Points 
on one side of the boundary represent perturbations which have a positive 
growth rate cr k , while points on the other side represent perturbations with 
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negative growth rates which are therefore suppressed in the solution. Thus 
we expect that solutions of the full differential equation will consist only of 
modes whose wave numbers are on the unstable side of the boundary. The 
onset of instability of the flat sandbed occurs at the value J c of J for which 
only a single mode, with wave number k c , is marginally stable and no other 
modes are unstable. These critical values may be determined by solving 
a = and da/dk = simultaneously, which yields 



Eliminating J c gives 



Jc 


sinfc c 


= 




'k c ., 


Jc 


cosk c 


= 




-1. 




ta.nk c 


— 


k c . 





(15) 



(16) 



Thus the critical values are computed as k c = 4.493 and J c = 4.603. The 
wavelength of the marginal mode (in units of Lq) is A c = 2ir jk c = 1.398, 
somewhat longer than the flight distance of a grain in saltation. For k = k c , 
the NO saltation relation leads to u c = —k c , so the phase velocity of the 
marginal mode is v = w c /k c = — 1. With the symmetric saltation relation we 
get v = — (1 + J c ) = —5.60. This is a surprising result of the model, that while 
the sand grains that form the ripples are blown downwind, the ripple pattern 
itself drifts upwind. The group velocity, however, is large and positive: From 
(ITJP we get duik/dk = J(—ksmk + cosk), which goes to k 2 — 1 = 19.19 at 
the critical point. For the symmetric saltation relation, the group velocity is 
lower by J, so at critical it is 14.58. Note that all velocities are in units of 
D/L . 

If we make the problem two-dimensional, allowing the sandbed to extend 
in both x and y directions, very little changes. The creep term in the evolu- 
tion equation becomes D V 2 /i, and as a result the expression for o^ changes 
to 

a(k,k y ) = —k 2 — ky — J ksink, (17) 

where k is now the x component of the wave vector of the Fourier mode and 
k y is its y component. Clearly, the linear growth rate for a mode with nonzero 
k y is always less than the rate for the corresponding mode with k y = 0. Thus 
we do not expect to see instabilities in which the transverse shape of the 
ripples becomes wavy, since the first instability to occur is against a mode in 
which the ripples are parallel to the y axis. 
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4 Nonlinear Analysis 

We now carry out a weakly nonlinear analysis to determine the amplitude, 
shape, and propagation velocity of the restabilized ripple patterns which form 
when J is slightly above its critical value J c . The nonlocality of the model, 
the dispersion in the imaginary part of the linear growth rate, and the lack 
of an up-down symmetry lead to some unusual features in the analysis. 

We begin with the assumption that the fundamental wave number k of 
the pattern which develops does not deviate much from the critical value k c 
when J is near J c . Hence we define a small parameter e by setting 

J = J c + e 2 , (18) 

and then define a scaled wave number deviation q by writing 

k = k c + eq. (19) 

This is the appropriate scaling for the wave number because the stability 
boundary is approximately quadratic in k — k c and linear in J near its max- 
imum. Substituting these expressions into the linear growth rate (|13"D and 
expanding to second order in e gives 

° k = \t 2 kl(j c -q 2 )+0{e\ (20) 

Furthermore, from the expression (jTJj) for uj we find that the phase velocity 
of the ripples is given by 

v = w /k = -l + ek c q - l -e 2 (j - g 2 ) + 0(e 3 ). (21) 

The first stage of the nonlinear analysis consists of expanding the evo- 
lution equation (|7|) in powers of h, assuming the overall amplitude of h is 
small. To do this, we may rewrite the NO saltation relation (|6]) in the form 

Z = x-l-h(£,t), (22) 

repeatedly substitute this expression for £ back into the h(£,t) on the right 
side, and finally expand in powers of h. Differentiating the result with respect 
to x then gives 

|L = 1 _ h '(x - 1, t) + \\h\x - 1, t)\" - l -[h\x - 1, £)]'" + 0{h'). (23) 



To third order in h, then, the evolution equation becomes 

^^ = h'\x,t)-jh'(x-i,t)+^-[h 2 (x-i,t)r~^[h 3 (x-i,t)r+... (24) 

Ot 2 o 

This is the equation whose ripple solutions we will presently compute. 

It is remarkable that the expansion of d^/dx has such an economical form. 
In fact it is not difficult to show that the pattern continues to all orders in 
h. To see this, consider the integral 

If= f°° f(x-l)^dx, (25) 

J — oo CtX 

where £(x) is given by the saltation relation (|6]) and / is a test function which 
is integrable and infinitely differentiable, but otherwise arbitrary. We now 
change variables in this integral from x to £, 

/oo 
f{Z + h{£))dZ, (26) 

-oo 

and expand the integrand in powers of h to get 

WgsF^"®* (27) 

Next we integrate the fcth term by parts k times to get 

and finally change variables again from £ to x = C, + 1 , 

This result has the same form as the original expression for If, but with d^/dx 
replaced by an expansion. However, since the test function / is arbitrary, 
this requires the expansion and d^/dx to be equal: 

dj _ f, (~l) fc d k h\x - 1) 

dx f^ k\ dx k ' l ' 

fe=0 



We now turn to the second stage of the calculation, namely finding so- 
lutions to the third-order approximation (]24f) to the evolution equation. We 
assume the solution will have a fundamental wave number k in the unstable 
range, with an amplitude of order e. The quadratic terms in the evolution 
equation will then generate a Fourier component in the solution with wave 
number 2k and possibly a constant term, and the cubic terms will lead to a 
component with wave number 3k. Thus we write 

h(x, t) « eM(t) cos(kx - <p(t)) + e 2 M (t) + e 2 M 2 (t) cos[2(fcr - <p{t)) + 6 2 {t)} 

+ e 3 M 3 (t) cos[3(A;x - 0(t)) + 6 3 (t)} + • • • ,(31) 

allowing phase differences among the various Fourier components. We sub- 
stitute this ansatz into the evolution equation and expand in powers of e. 
Then the coefficients of cos(fcr — </>(£)) and sin(foc — </>(£)) give equations for 
the fundamental amplitude M(t) and phase <p(t): 

M = a k M - e 2 Jk 2 M M cos k Jk 2 M 2 Mcos(k - 2 ) 

+-Jk 3 M 3 smk + 0(e 4 ) (32) 

8 

2 2 

= Ci ; fc - e 2 JA; 2 MosinA;--JA; 2 M2sin(A;-0 2 )--JA; 3 M 2 cosfc + O(e 4 ). (33) 

2 8 

Evidently we need to find Mo, M 2 , and #2 in order to determine the am- 
plitude M and propagation velocity <p/k of the pattern. The x-independent 
term in the expansion of the evolution equation gives M = 0, so M Q is in 
fact a constant; as argued above, we can choose it to be zero, so the M 
terms in the M and <fi equations can be dropped. The equations for M 2 and 
62 come from the coefficients of cos(2fcx — 2^ + 62) and sin(2fcx — 2^ + 62) in 
the evolution equation. These are best written in the form 

^M 2 e-* = (a 2k + iu 2k - 2too k )M 2 e-^ - Jk 2 M 2 e 2ik + 0(e 2 ). (34) 
at 

Similarly, we find equations for M 3 and 8 3 , 

—M 3 e- i0i = (a 3k +iu 3k -3iu k )M 3 e- Wi --Jk 2 MM 2 e 3ik - ie2 --iJk 3 M 3 e 3ik +0(e 2 ) 
at 2 8 

(35) 
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Note that in the equation for M(t), all the terms on the right are of order 
e 2 . Therefore M(t) changes on a long time scale of order e~ 2 , while M 2 and 
6 2 vary on times of order unity. Thus we may regard M as a constant in 
the equations for M 2 and M 3 . Since a 2 k is negative, M 2 exp(— i6 2 ) goes to a 
quasi-steady state value which is proportional to M 2 . Substituting this value 
into the M equation gives 

M = a k M - e 2 AM 3 + 0(e 4 ) (36) 

The analytical expression for the Landau constant A is complicated and un- 
enlightening; substituting (|18|) and (|19|) into it gives 

A = 16.905 + 64.680eg + O(e 2 ) (37) 

Note that the correction terms in the evolution equation for M, which would 
come from including higher-order terms in the original expansions (E3j) for 
the evolution equation and (|31~D for h(x,t), are of order e 4 , not e 3 . As a 
result, the e 3 terms in the equation come only from expanding the analytical 
expressions for a and A in powers of e. This also holds for the equations for 
the higher harmonics. Thus we get the first-order corrections to all of our 
results essentially for free. 

From (37) and the third-order expansion ( P0"D for 0^. we obtain the steady- 
state amplitude M to first order in e, 

M 2 = (0.25946 - 0.59719g 2 ) - (0.87725 - 2.10774g 2 ) tq + 0(e 2 ) (38) 

We then find the phase velocity, 

v ph = tu/k = -l+4.4934eg+(1.877-4.320g 2 )e 2 -(3.439-10.128g 2 )e 3 g+O(e 4 ), 

(39) 
and the group velocity, 

v gr = du/dk = 19.1907+13.4802eg+(18.241-40.984g 2 )e 2 -(43.198-107.23g 2 )e 3 g+O(e 4 

(40) 
from (33). The equation (34) for M 2 and 9 2 gives 

M 2 /M 2 = 0.908115 + 0.50118eg, 

9 2 /n = 0.22916 - 0.36189eg, (41) 
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and from (35) we find 

M 3 /M 3 = 1.52273 + 1.89385eg, 

9 3 /n = 0.4309 - 0.6024eg (42) 

As usual, the ripple solutions we have found are not all stable; instead, 
those with too large a wave number deviation q are linearly unstable. The 
calculation of the critical value of q is rather intricate, so we defer it to the 
Appendix. The result is that the range of stable wave numbers is rather 
wider than usual - it extends out to q — 0.9095g , where q = (2/ J c ) 1 ^ 2 is 
the wave number deviation at which a k vanishes to leading order in e. At 
the edge of the stable range, the amplitude M of the ripple solution is 0.4157 
times its value ay q = 0. 

If instead of the NO saltation relation we use the symmetric relation (§), 
the results of the analysis are rather different. The expansion of d^/dx is not 
as simple and clean as the derivation above; the evolution equation (|2~4D is 
replaced by 

a/t fo*) = h"{x,t)-J[h{x-l,t)-h{x,t)]' + J{[h{x-l,t)-h{x,t)]ti{x-l,t)y 

-{[h(x-l, t)-h{x, t)] 2 h"(x-l, t)}'-J{[h{x-l, t)-h(x, t)} [ti{x-l, t)] 2 }'+- ■ ■ 

(43) 
We again substitute the ansatz fl3"T| ) into this equation and work out the 
Fourier components of the result. The equations for M and </> become 

M = a k M - e 2 Jk 2 M 2 M[cos k cos 9 2 - cos(2A; - 2 )] 

+ — Jk 3 M 3 (sin k -2 sin 2k), (44) 

e 2 

(fy = u k + e 2 Jk 2 M 2 [cos k sin 6 2 + sin(2k- 6 2 )] Jk 3 M 2 (cos k- cos 2k) (45) 

where now uo k is given by —Jk (1 — cos A;) as is appropriate for this model. 
Note that the Mq terms which were present in (32) and (33) are absent 
here; this is because the new saltation relation respects the symmetry under 
addition of a constant to h. The evolution of M 2 and 62 is now given by 

j t M 2 e- W2 = [a 2k + iuj 2k - 2iu k )M2e-^ + Jk 2 {e ik - e 2ik )M 2 + 0(e 2 ), (46) 
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and the third harmonic by 

jMze-** = (a 3k + tcu 3k -3 l cu k )M 3 e-^ + pk 2 e^(l-e tk )(l + 3e tk )MM 2 e' l(>2 

+ li.Jk 3 e tk (l - e lk ) (1 - 3e tk )M 3 + 0(e 2 ) (47) 

8 

After carrying out the calculation we find a much larger value for the 

Landau coefficient, 

A = 151.26 + 88.014eg + O(e 2 ). (48) 

Thus for a given wave number, the restabilized amplitude M of the ripples 
is smaller by a factor of about 3: 

M 2 = (0.0289975 - 0.066743g 2 ) - (0.003966 - 0.0190315g 2 ) tq + 0(e 2 ). (49) 

The phase velocity is more negative than before, as we found from the linear 
stability analysis, 

v ph = -5.6033+4.4934eg-(1.506-1.16475g 2 )e 2 -(0.197-1.853g 2 )e 3 g+O(e 4 ), 

(50) 
and likewise the group velocity, 

Vgr = 14.5874+13.4802eg-(2.569-4.612g 2 )e 2 -(8.152-19.799g 2 )e 3 g + O(e 4 ). 

(51) 
For a given amplitude, however, the harmonics are stronger than before: we 
find 

M 2 /M 2 = 1.41691 + 0.213856eg, 

2 /7r = rrO.4443 - 0.2027eg, (52) 

and 

M 3 /M 3 = 2.89266 + 1.57014eg, 

3 /n = -0.1426 - 0.3952eg. (53) 

The range of wave numbers for which these solutions is stable is some- 
what narrower than before but still wider than usual, extending out to 
q = 0.7571g , where the amplitude M is 0.6533 times its value at q — 0. 
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5 Numerical solutions 

We now present numerical solutions and compare them with the predic- 
tions made in the previous section. The nonlocal evolution equation (0) 
was solved numerically with periodic boundary conditions on a system of 
length I = 2ir/k, so that only the Fourier modes nk contributed to the 
solutions. For the discretization scheme, we chose an explicit method us- 
ing forward differences in time and central differences in space. The axis 
was discretized at 2 9 equally spaced sites with Ax = 2ir/k • 2 -9 and solu- 
tions were generated for five different values of J near J c = 4.603, namely 
J = 4.62,4.65,4.70,4.75,4.80, and values of k were chosen to span the un- 
stable region. Initial conditions were sinusoids of wave number k centered 
around h = 0. At t=0, we start with a sinusoid at a particular wave number 
k, and let it evolve with (Ax) 2 / At = 1/4 until it reaches a steady state. 
This takes about 10 6 time steps. The nonlocal term in the evolution equa- 
tion, J (d^/dx — 1), was evaluated for a given x by finding the nearest upwind 
value of £ satisfying the equation x = £ + 1 + h(£). Spefifically, the first root 
of the function /(£; x) = £ + 1 + /i(£) — x with value less than x was obtained 
by simply finding the two sites upwind of x and nearest to it between which 
/(£; x) changed sign. Then, d^/dx was calculated using the values of h at 
these sites. The final steady state, h(x,t), is then Fourier transformed, i.e.: 

oo 

h(x,t) = y^[q n (t) smnkx + b n (t) cos nkx], (54) 

n=l 

from which we obtain M n = (o^ + b 2 /) 1 ^ 2 and 9 n = tan _1 (6 n /a n ). Note that 
Mi = M and 9\ = <fi in this notation. The nonlinear analysis in the previous 
section predicts that these quantities will go to time-independent values. We 
find numerically that they actually oscillate as a function of time around 
their mean values. However, the magnitudes of these oscillations are quite 
small and decrease with increasing grid resolution, so we believe them to be 
numerical artifacts. We therefore take the time averaged mean values and 
compare them with the predictions of the weakly nonlinear analysis. 

We also find that although we start with an initial profile with the average 
height h = 0, the mean position of the steady state pattern shifts slightly 
upward in some cases, downward in others, to a small but finite ho- Since 
the mean height of the sandbed is conserved by the exact evolution equation, 
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we believe that this is also a numerical artifact. Moreover, as mentioned 
in section 2, we can map any steady state solution with finite ho to the 
solution with /i = by redefining the horizontal length scale from L to 
L (l + ho) and shifting the control parameter from J to J(l + ho). However, 
the magnitude offset ho was always of the order of 10~ 5 to 10~ 3 , and thus in all 
cases studied here, the corrections due to such an offset are quite negligible. 
Hence, our results without these corrections are virtually identical to those 
with corrections. 

Figure 2 shows the amplitude M of the fundamental mode as a function 
of k for different values of J. The data points are fairly close to the values 
predicted by the first-order expansion (38), which are represented by contin- 
uous curves. Note that the curves are asymmetric around the critical value 
k c and the asymmetry becomes more pronounced for larger J. The weakly 
nonlinear analysis is capable of predicting this asymmetry only because the 
order-e terms are included. 

In Figure 3 we plot the phase velocity of the fundamental mode against 
k for different values of J . The speed was obtained by calculating <f){t) in the 
expression cos(kx — <p{t)), which is proportional to the fundamental mode 
in the steady state. The function <p(t) was found to be linear in t, so v was 
calculated as v — d<f)(t)/dt/k. The data points are compared against the 
weakly nonlinear predictions (solid line) given by (39). Only for fairly large 
J, and only near the high-wave-number end of the band of ripple solutions, 
does the velocity become positive, that is, in the direction of the wind. 

In Figures 4a and 4b are plotted the ratios of the amplitudes of the 
second and third harmonics to the appropriate powers of the fundamental 
ampltude, i.e., R 2 = M 2 /M 2 and R 3 = M 3 /M 3 . The data fit quite well 
with the theoretical predictions for both cases, in particular near the onset 
k c = 4.493, where the nonlinear analysis is most reliable. Note that the 
first-order terms in the analytical results match the slope of the numerical 
results. The curvature which is evident in the numerical data for k farther 
from k c is apparently a higher-order effect. Note that the width of the band 
of ripple solutions increases with e, so an appreciably large e is required to 
reach these larger values of \k — k c \. 

In Figures 5a and 5b we plot the phase angles 6 2 and #3 against k. The 
agreement between the simulations and the weakly nonlinear analysis is again 
quite strong for J near onset. The order-e terms in the analytical results 
match the slope of the numerical data. For higher J we observe a systematic 
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downward deviation in the numerical results. The shift appears to be linear 
in J, and so is second-order in e = (J — J c ) 1 ^ 2 - 

6 Discussion 

We have carried out numerical and weakly nonlinear analyses of the Nishimori- 
Ouchi continuum model [3,11] for windblown sand, and also for a modifica- 
tion of that model which respects the physical symmetry of the system under 
changes of the reference level of height. Both versions of the model yield the 
surprising result that the ripple patterns, which form when the flat sandbed 
becomes unstable, drift upwind even as the sand which forms the ripples 
is blown downwind. This drift is found in the linear stability analysis and 
persists in the weakly nonlinear results, and numerical integrations confirm 
that it is a real consequence of the model. Such a counterintuitive result 
has not been examined or detected by previous Monte Carlo simulations of 
this model [3,11] or in real experiments [10]. It would be interesting to check 
experimentally whether or not ripples can move against the wind. The sym- 
metric version of the model actually predicts a considerably higher upwind 
drift speed than the original Nishimori-Ouchi version. 

It may also be surprising that the differences between the symmetric 
and Nishimori-Ouchi models are merely qualitative. The restabilized ripple 
pattern for a given value of the control parameter have smaller amplitudes 
(by a factor of about 3) and higher drift velocities (by a factor of over 5) 
in the symmetric model than in the original version. The relative sizes and 
phases of the higher harmonics in the ripple shape are also different for the 
two models. 

A number of modifications to the model are needed in order to make 
it comparable with experiments. A major ingredient that is left out of the 
model is any effect of the surface topography on the wind. This lack means 
that there is no shadowing effect in the model. Including such an effect would 
make it more likely for grains to settle on the downwind side of a ripple than 
on the upwind side, and more likely for them to be blown off the upwind 
side than the downwind side. This would likely reduce the tendency of the 
ripples to drift upwind. The result that the ripples drift upwind in this model, 
which neglects shadowing, may be an indirect indication of the importance 
of shadowing in the development of real ripple patterns. An improved model 
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of creep may also be needed; a downwind bias in the creep would modify the 
drift velocity. Perhaps most critical is a better and more realistic form of the 
saltation function, which must account the effects of the topography of the 
sandbed, and the many particle dynamics of the grains in the air as well as 
on the surface. 
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A Stability of ripple solutions and the Eck- 
haus boundary 

In this section, we examine the stability of the ripple solutions, and so deter- 
mine the Eckhaus boundary in the J-k plane, within which the solutions are 
linearly stable and so may be observed. We begin with the ripple solution 

ho(x,t) = eM cos(kx - ut) + e 2 M 2 cos[2(kx - ut) + 6 2 ] + 0(e 3 M 3 ), (A- 1) 

with k = k c + eg, and add an infinitesimal perturbation hi(x,t). If the 
perturbation contains a Fourier component with wave number k + eg', then 
the nonlinear terms in the evolution equation will generate a component with 
wave number k — eg'. Thus we will start the calculation by taking h\ to have 
the form 

hi(x, t) = A_(t) cos[(k-eq')x-(u)t+(f)-(t))]+A + (t) cos[(k+eq')x-(u)t+(j) + (t))]. 

(A-2) 
Substituting this into the evolution equation, expanding, and picking off the 
coefficients of the sines and cosines of (k — eq')x and (k + eq')x yields a closed 
set of equations for the amplitudes A_ and A + and the phase + + 0_. These 
equations have the form 

A_ = £_A_ + aA + cosijj, 

A + = Y* + A + + aA^ cosip, 



ij> = n- a[{A + /A_) + {A_/A + )\ sin^ 



(A -3) 



where ip is + + 0_ plus a constant which depends on k and q' (but not time), 
and the overdot denotes a derivative with respect to the slow time variable 
e 2 t. The coefficients are given to leading order in e by 



k 2 



T-f) M2 ^' 



Q' 



a 



a k\ M 2 
1 

2 



n 



1.98183 k 2 c M 2 , 
k 3 c M 2 + 3k c q' 2 . 



{A -A) 



19 



Note that it is important to keep the second-order term in ho during the 
calculation, since it contributes to a. (Omitting it changes ceo to 2.58559, a 
30% change.) 

We must now determine whether the amplitudes given by (A-3) grow or 
decay with time. We can simplify the equations somewhat by defining 

R = A + /A^; (A - 5) 

the amplitude equations then become 

A_ = (E_ + aRcosip)A + , 

R=(H+- E_)i? + a(l - R 2 ) cos4>, 

ip = Q- a[(l + R 2 )/R] sin i). (^-6) 

Note that the first equation decouples from the last two. If it happens that 
R and ip go to constants as t — *■ oo, then the amplitudes decay for £_ + 
aRcosip < and grow otherwise. Thus for a given q, the ripple state (A- 
1) is linearly stable if this inequality is satisfied for every q', otherwise it 
is unstable. To see what R and ip actually do, we combine the R and ijj 
equations into an evolution equation for the complex variable 

Z = Rexp(i^j), (A -7) 

namely 

Z = «(l-Z 2 ) + (S + -S_+ifi)Z. (-4-8) 

Clearly this equation has two fixed points, and solving it exactly reveals that 
the one with positive real part is a global attractor and the one with negative 
real part a global repeller. Thus Rcosip does go to a constant, and from its 
value we can decide whether the ripple state is stable or not. 

Since it is the real part of Z, namely Rcosip, which determines whether 
the perturbation grows or decays, it is useful to rewrite (A-8) in terms of the 
real and imaginary parts of Z, 

Z = X + iY. (A -9) 

We find 

X = a(l - X 2 ) + (E + - S_)X - ttY + aY 2 , 
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Y = (S+-S_)r + fiX-2aXF. (^-10) 

By linearizing about a fixed point (X, Y) of this system, we quickly find that 
the fixed point is an attractor for X > (E + — E_)/2a. To find the fixed 
points, we set X = Y = and solve the second equation for Y in terms of 
X, then substitute into the first equation to get 

The two sides of this equation are plotted in Fig. 6. Both sides are symmetric 
about X = (E + — XL) /2a, so there is clearly one solution with X greater 
than this - the attractor - and one, the repeller, with X less. In order for the 
perturbation to decay, the attractor must have X less than — E_/a. From 
the plot, we see that this means that at X = — E_/a the right side of (A-ll) 
must be greater than the left side. After a little algebra, we can write this 
condition for the perturbation to decay in the form 

E E , a 2 (E + + E.)3 

Equation (A- 12) above is the condition for the amplitude of a perturba- 
tion with a specific value of q' to decay. In order to conclude that the ripple 
solution with a given q is linearly stable, we must see to it that this condition 
is satisfied for all q'. For this we must substitute for the parameters from 
(A-4) above. To put the result into a useful form, we define Q = 2q' 2 /k 2 M 2 
and eliminate q 2 in favor of M 2 . After some rearranging, we find that the 
condition for the solution to be stable is 

M 2 16Q k 2 (P + Q) 2 + (1 + 3Q) 2 

J k 2 c [{(3 - Q) 2 + AQ] [k 2 c {(3 + Q) 2 + (1 + 3Q) 2 ] - \Qa 2 {(3 + Q) 2 ' 

(A- 13) 
where (3=1 — (4A/A; 2 ) = 0.834132. The complicated function of Q on the 
right has a single maximum for positive Q, at a height of 0.04484. Ripple 
states with M 2 below this are unstable, while those with M 2 larger than 
this are linearly stable. From this we find that the range of wave numbers of 
linearly stable ripple solutions is given by \q\ < 0.9095go, where go = (2/ J c ) 1 ^ 2 
is the largest wave number for which a ripple solution exists. 
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In summary, we have found that in the weakly nonlinear regime, the flat 
sandbed is unstable against perturbations with wave numbers k in the range 



| /,•/,-,!< a /0 = \/2( j -I) U- 11) 



while the Eckahus boundary is given by 



\k - k c \ < 0.9095eg = Jl.654 (— - l). {A- 15) 

For the symmetric saltation relation, the structure of the calculation is 
the same but the numbers are different. We find that the marginally stable 
wave number is given by q = 0.7571go, so the Eckhaus boundary is now given 
by 

\k - fe c | < 0.7571eg = W 1.1465 (— - l\ {A- 16) 
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Figure Captions 



Fig. 1: Saltation refers to the process of a single grain being ejected from 
the surface at a point £ and being blown to a landing point x by the wind. 
Fig. 2: The steady state amplitude eM vs. k for five different values of J . 
The continuous lines are the analytical predictions given by Eq. (38). 
Fig. 3: The propagation speed of the steady state patterns v p h vs. k for 
different values of J. Numerically obtained values are compared to the ana- 
lytical predictions (continuous curves) from Eq. (39). 

Fig. 4: The ratios (a) R 2 = M 2 /M 2 and (b) R 3 = M 3 /M 3 are plotted against 
k. The solid line is the analytical predictions from Equations (41) and (42). 
Fig. 5: (a) 9 2 and (b) 83 are plotted against k for five different values of J. 
The solid line is the analytical predictions from Equations (41) and (42). 
Fig. 6: The right and left hand sides of Eq. (A-ll). 
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Typical Saltation Trajectory and Related Variables 
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